SOFTWARE PARA EL ANÁLISIS DE DATOS (SAD)

MÁSTER UNIVERSITARIO EN BIOINFORMÁTICA Y BIOESTADÍSTICA

Motivación y datasets empleados

El objetivo de nuestro trabajo es estudiar si existe alguna relación entre la vacuna BCG (Bacillus de Calmette y Guérin) para la tuberculosis y los datos de mortalidad de la COVID-19 en algunos países, ya que hay estudios que sugieren esta vacuna incrementa las capacidades inmunitarias de la población, hecho que se ve en el número reducido de fallecimientos por COVID-19 en ciertos países. Mediante los conjuntos de datos de BCG y de mortalidad por COVID-19 cedidos por The BCG world atlas y por BCG - COVID-19 AI Challenge de Kaggle, vamos a intentar desvelar dichas relaciones.

Los ficheros en cuestión son del tipo `csv, así que son fácilmente importables a data frames enR:

# Cargamos ambos datasets,

BCG_strain <-
  read_csv("task_2-BCG_strain_per_country-1Nov2020.csv")

COVID_noformat <-
  read_csv(
    "task_2-COVID-19-death_cases_per_country_after_fifth_death-till_22_September_2020.csv"
  )

# Intenté ver que hay dentro de los data frames, pero el print es feo así que lo escribiré a mano
# str(COVID_noformat)
# str(BCG_strain)

El contenido de las variables BCG_strain y COVID_noformat es el siguiente:

BCG_strain COVID_noformat
country_name country_name
country_code alpha_3_code
mandatory_bcg_strain_2015-2020 date_first_death
mandatory_bcg_strain_2010-2015 date_fifth_death
mandatory_bcg_strain_2005-2010 deaths_per_million_10_days_after_fifth_death
mandatory_bcg_strain_2000-2005 deaths_per_million_15_days_after_fifth_death
mandatory_bcg_strain_1990-2000 deaths_per_million_20_days_after_fifth_death
mandatory_bcg_strain_1980-1990 deaths_per_million_25_days_after_fifth_death
mandatory_bcg_strain_1970-1980 deaths_per_million_30_days_after_fifth_death
mandatory_bcg_strain_1960-1970 deaths_per_million_35_days_after_fifth_death
mandatory_bcg_strain_1950-1960 deaths_per_million_40_days_after_fifth_death
vaccination_timing_unified deaths_per_million_45_days_after_fifth_death
BCG Atlas: Which year was vaccination introduced? deaths_per_million_50_days_after_fifth_death
Year of changes to BCG schedule deaths_per_million_55_days_after_fifth_death
BCG Atlas: BCG Recommendation Type deaths_per_million_60_days_after_fifth_death
BCG Atlas: Details of changes deaths_per_million_65_days_after_fifth_death
BCG Atlas: Timing of 1st BCG? deaths_per_million_70_days_after_fifth_death
BCG Atlas: BCG Strain deaths_per_million_75_days_after_fifth_death
BCG Atlas: How long has this BCG vaccine strain been used? deaths_per_million_80_days_after_fifth_death
deaths_per_million_85_days_after_fifth_death
deaths_per_million_90_days_after_fifth_death
deaths_per_million_95_days_after_fifth_death
deaths_per_million_100_days_after_fifth_death
deaths_per_million_105_days_after_fifth_death
deaths_per_million_110_days_after_fifth_death
deaths_per_million_115_days_after_fifth_death
deaths_per_million_120_days_after_fifth_death
deaths_per_million_125_days_after_fifth_death
deaths_per_million_130_days_after_fifth_death
deaths_per_million_135_days_after_fifth_death
deaths_per_million_140_days_after_fifth_death
deaths_per_million_145_days_after_fifth_death
deaths_per_million_150_days_after_fifth_death
stringency_index_10_days_after_fifth_death
stringency_index_15_days_after_fifth_death
stringency_index_20_days_after_fifth_death
stringency_index_25_days_after_fifth_death
stringency_index_30_days_after_fifth_death
stringency_index_35_days_after_fifth_death
stringency_index_40_days_after_fifth_death
stringency_index_45_days_after_fifth_death
stringency_index_50_days_after_fifth_death
stringency_index_55_days_after_fifth_death
stringency_index_60_days_after_fifth_death
stringency_index_65_days_after_fifth_death
stringency_index_70_days_after_fifth_death
stringency_index_75_days_after_fifth_death
stringency_index_80_days_after_fifth_death
stringency_index_85_days_after_fifth_death
stringency_index_90_days_after_fifth_death
stringency_index_95_days_after_fifth_death
stringency_index_100_days_after_fifth_death
stringency_index_105_days_after_fifth_death
stringency_index_110_days_after_fifth_death
stringency_index_115_days_after_fifth_death
stringency_index_120_days_after_fifth_death
stringency_index_125_days_after_fifth_death
stringency_index_130_days_after_fifth_death
stringency_index_135_days_after_fifth_death
stringency_index_140_days_after_fifth_death
stringency_index_145_days_after_fifth_death
stringency_index_150_days_after_fifth_death

Una visualización preliminar de estos datos revela que son todos del tipo string y que además muchas columnas sin datos (columnas cuyo único contenido es NULL), por lo tanto llevaremos a cabo una limpieza de los mismos además de cambios de tipo de variables para que las manipulaciones posteriores sean más cómodas. Los detalles se muestran en el siguiente bloque de código:

# Limpiar datos de BCG

# Elimino columnas que sean sólo NA
BCG_strain <- BCG_strain[, apply(!is.na(BCG_strain), 2, all)]

# De momento, no me interesa qué vacunas se ponían cada año, sino si se ponían o no.

# Transformo los valores de cada año en
# 0 - No se ponía vacuna, hasta ahora None
# 1 - Sí se ponía vacuna
# NA - Este dato es desconocido, hasta ahora Unknown

BCG_strain_no_strain <- BCG_strain
# Transformo los valores de las columnas
BCG_strain_no_strain[, -1] <-
    sapply(BCG_strain_no_strain[, -1], function(x) {
        a <-
            gsub("None", 0, x) %>% gsub("Unknown", NA, .) # Añado los 0 y los NA.
        for (i in 1:length(a)) {
            # Serán 1 aquellos que no sean ni 0 ni NA
            if (a[i] != "0" && !is.na(a[i])) {
                a[i] <- 1
            }
        }
        return(as.integer(a)) # Cambio las columnas a integer
    })
BCG_no_strain_no_NA <- na.omit(BCG_strain_no_strain)

#Versiñon más compacta del dataframe, sin datos a diferentes días o años.
#Agrupando los datos de vacunas en tres columnas:
#periods_with_vaccine - incluye el número de periodos estudiados con vacunación activa
#vaccination_2020_2015 - el único periodo con el que nos quedamos, el último
#first_vaccine_year - de los años estudiados, el primero con campaña de vacunación. En el caso de no tener vacunación, este será el último año estudiado (2020)
#last_vaccine_year - de los años estudiados, el último con campaña de vacunación. En el caso de no tener vacunación, este será el primer año estudiado (1950)

#Creamos el nuevo dataframe simplificado
BCG_no_strain_simple <- data.frame(
  "country_name" = BCG_no_strain_no_NA$country_name,
  "periods_with_vaccine" = BCG_no_strain_no_NA%>%
  .[2:ncol(.)] %>% rowSums(), #sumamos los periodos con vacuna
  "vaccination_2020_2015" = BCG_no_strain_no_NA$`mandatory_bcg_strain_2015-2020`)

#Añadimos el último año con vacunación
BCG_no_strain_simple$last_vaccine_year <- names(BCG_no_strain_no_NA[2:ncol(BCG_no_strain_no_NA)])[max.col(BCG_no_strain_no_NA[2:ncol(BCG_no_strain_no_NA)] != 0, 'first')] %>%
  substring(nchar(.)-3, nchar(.))%>%
  as.numeric()

#Añadimos el primer año con vacunación
BCG_no_strain_simple$first_vaccine_year <- names(BCG_no_strain_no_NA[2:ncol(BCG_no_strain_no_NA)])[max.col(BCG_no_strain_no_NA[2:ncol(BCG_no_strain_no_NA)] != 0, 'last')] %>%
  substring(nchar(.)-8, nchar(.)-5)%>%
  as.numeric()


#El próximo código es necesario para que los países sin campaña de vacunación no obtengan los mejores valores. El código utilizado para obtener el último o primer año de vacunación les favorece, ya que obtiene el primer o el último índice de aquellos valores distintos de 0. Como en su caso no hay ningún valor distinto a 0, este sería simplemente el primero o el último. Para arreglar esto, establezco manualmente que tengan el último año de vacunación más bajo posible y el primer año de vacunación más alto posible.
BCG_no_strain_simple[BCG_no_strain_simple$periods_with_vaccine == 0,]$last_vaccine_year = 1950

BCG_no_strain_simple[BCG_no_strain_simple$periods_with_vaccine == 0,]$first_vaccine_year = 2020
####################################################################################

# Limpiar datos de COVID

# Elimino columnas que sean sólo NA
COVID_noNA <- COVID_noformat[, apply(!is.na(COVID_noformat), 2, all)]

# En este caso, para variar, los valores vacíos están denotados como NULL,
# cambio esto a NA
COVID_Na <- sapply(COVID_noNA, function(x)
    gsub("NULL", NA, x))

# El resulatado de la función anterior es una string. Lo convierto a dataframe.
COVID_Na_df <- as.data.frame(COVID_Na)

# Modifico las fechas para que se almacenen como Date
COVID_Na_df[, c("date_fifth_death")] <-
    as.Date(COVID_Na_df[, c("date_fifth_death")], "%d/%m/%y")

COVID_Na_df[, c("date_first_death")] <-
    as.Date(COVID_Na_df[, c("date_first_death")], "%d/%m/%y")

# Modifico las muertes para que se almacenen como floats.
COVID_Na_df[, -c(1, 2, 3, 4)] <-
    sapply(COVID_Na_df[, -c(1, 2, 3, 4)], as.numeric)


COVID_Na_df2 <- data.frame("country_name" = COVID_Na_df$country_name,
                           "dpm_100d" = COVID_Na_df$deaths_per_million_100_days_after_fifth_death,
                           "si_100d" = COVID_Na_df$stringency_index_100_days_after_fifth_death)
################################################################################################

# Junto ambos dataframes en uno sólo.
COVID_BGC <-
    inner_join(BCG_no_strain_simple, COVID_Na_df2, by = "country_name")

Nuestra tabla resultante es la siguiente:

Tabla 1. Vacunación de BCG por países y muertes por COVID-19
country_name periods_with_vaccine vaccination_2020_2015 last_vaccine_year first_vaccine_year dpm_100d si_100d
Afghanistan 6 1 2020 1980 26.00 78.70
Angola 6 1 2020 1980 4.38 NA
Argentina 3 0 2010 1990 30.64 92.59
Armenia 4 1 2020 2000 190.67 NA
Australia 9 1 2020 1950 4.00 52.31
Bangladesh 6 1 2020 1980 11.95 74.54

Mediante esta tabla llevaremos a cabo nuestros análisis. A continuación mostramos la estructura de la misma:

str(COVID_BGC)
## 'data.frame':    65 obs. of  7 variables:
##  $ country_name         : chr  "Afghanistan" "Angola" "Argentina" "Armenia" ...
##  $ periods_with_vaccine : num  6 6 3 4 9 6 6 9 8 9 ...
##  $ vaccination_2020_2015: int  1 1 0 1 1 1 1 1 1 1 ...
##  $ last_vaccine_year    : num  2020 2020 2010 2020 2020 2020 2020 2020 2020 2020 ...
##  $ first_vaccine_year   : num  1980 1980 1990 2000 1950 1980 1980 1950 1960 1950 ...
##  $ dpm_100d             : num  26 4.38 30.64 190.67 4 ...
##  $ si_100d              : num  78.7 NA 92.6 NA 52.3 ...

Hagamos unos análisis descriptivos:

a) ¿Cuál es la media para la variable “deaths per million 150 days after fifth death” (representado por `dpm_100d) para el conjunto de países que la han medido? ¿Y para la variable “stringency index 150 days after fifth death”?

m1 <- mean(na.omit(COVID_BGC$dpm_100d))
m2 <- mean(na.omit(COVID_BGC$si_100d))
cat(sprintf("La media para dpm_100d es %.2f\n", m1))
## La media para dpm_100d es 110.44
cat(sprintf("La media para si_100d es %.2f\n", m2))
## La media para si_100d es 57.84

b) ¿Cuáles son los países cuyo valor para “deaths per million 100 days after fifth death” (en caso de estar presente) se encuentra por debajo de la media? ¿Y para la variable “stringency index 100 days after fifth death”?

cat(
  "Los países cuyo valor de dpm_100d es menor que la media del dataset son:\n\n",
  paste(subset(COVID_BGC, dpm_100d < m1)$country_name, collapse = ', '),
  "\n\n"
)
## Los países cuyo valor de dpm_100d es menor que la media del dataset son:
## 
##  Afghanistan, Angola, Argentina, Australia, Bangladesh, Bosnia and Herzegovina, Bulgaria, Colombia, Czech Republic, Denmark, El Salvador, Estonia, Ethiopia, Finland, Greece, Guam, Hungary, India, Indonesia, Japan, Kazakhstan, Kenya, Kuwait, Latvia, Malaysia, Malta, Nepal, Nigeria, Pakistan, Philippines, Poland, Saudi Arabia, Senegal, Sierra Leone, South Africa, Sudan, Taiwan, Tanzania, Thailand, Tunisia, Ukraine, Uruguay
cat(
  "Los países cuyo valor de si_100d es menor que la media del dataset son:\n\n",
  paste(subset(COVID_BGC, si_100d < m2)$country_name, collapse = ', ')
)
## Los países cuyo valor de si_100d es menor que la media del dataset son:
## 
##  Australia, Bosnia and Herzegovina, Bulgaria, Czech Republic, Estonia, Finland, Greece, Hungary, Ireland, Italy, Japan, Latvia, Malaysia, Poland, Senegal, Sierra Leone, Spain, Sweden, Switzerland, Taiwan, Tanzania, Thailand, Tunisia, Ukraine, Uruguay

c) ¿Cuáles son los países que cumplen ambas condiciones del apartado anterior?

cat(
  "Los países cuyos valores de dpm_100d y de si_100d son menores que la media del dataset son:\n\n",
  paste(subset(COVID_BGC, (dpm_100d < m1) &
                 (si_100d < m2))$country_name, collapse = ', ')
)
## Los países cuyos valores de dpm_100d y de si_100d son menores que la media del dataset son:
## 
##  Australia, Bosnia and Herzegovina, Bulgaria, Czech Republic, Estonia, Finland, Greece, Hungary, Japan, Latvia, Malaysia, Poland, Senegal, Sierra Leone, Taiwan, Tanzania, Thailand, Tunisia, Ukraine, Uruguay

d) ¿Cuáles son los países que han tenido campaña de vacunación de la vacuna BCG más reciente y que su la media de mortalidad a los 150 días es menor que la media?

cat(
  "Los países cuyos valores de dpm_100d son menores que la media del dataset y que además han tenido una reciente campaña de vacunación de la vacuna BCG son:\n\n",
  paste(subset(
    COVID_BGC, (dpm_100d < m1) &
      (vaccination_2020_2015 == 1)
  )$country_name, collapse = ', ')
)
## Los países cuyos valores de dpm_100d son menores que la media del dataset y que además han tenido una reciente campaña de vacunación de la vacuna BCG son:
## 
##  Afghanistan, Angola, Australia, Bangladesh, Bosnia and Herzegovina, Bulgaria, Colombia, El Salvador, Estonia, Ethiopia, Hungary, India, Indonesia, Japan, Kazakhstan, Kenya, Kuwait, Latvia, Malaysia, Malta, Nepal, Nigeria, Pakistan, Philippines, Poland, Saudi Arabia, Senegal, Sierra Leone, South Africa, Sudan, Taiwan, Tanzania, Thailand, Tunisia, Ukraine

Hagamos un poco de estadística descriptiva:

a) Resumen estadístico de las variables

# Obvimente hay variables en las que no tiene sentido hacer resumen estadístico, como el alpha_3_code, las strains... Pero por ahora lo voy a dejar
summary(COVID_BGC)
##  country_name       periods_with_vaccine vaccination_2020_2015
##  Length:65          Min.   :0.000        Min.   :0.0000       
##  Class :character   1st Qu.:5.000        1st Qu.:1.0000       
##  Mode  :character   Median :6.000        Median :1.0000       
##                     Mean   :6.277        Mean   :0.7538       
##                     3rd Qu.:8.000        3rd Qu.:1.0000       
##                     Max.   :9.000        Max.   :1.0000       
##                                                               
##  last_vaccine_year first_vaccine_year    dpm_100d          si_100d     
##  Min.   :1950      Min.   :1950       Min.   :  0.294   Min.   :19.44  
##  1st Qu.:2020      1st Qu.:1950       1st Qu.:  9.524   1st Qu.:38.89  
##  Median :2020      Median :1960       Median : 30.134   Median :59.72  
##  Mean   :2012      Mean   :1969       Mean   :110.444   Mean   :57.84  
##  3rd Qu.:2020      3rd Qu.:1980       3rd Qu.:116.554   3rd Qu.:74.07  
##  Max.   :2020      Max.   :2020       Max.   :584.686   Max.   :92.59  
##                                       NA's   :9         NA's   :12
# Igual habría que dibujar histogramas de tasas de muerte, y de vacunas puestas. No sé la verdad.

b) Distribución de las variables

Observemos cómo se distribuyen los datos de “deaths per million 100 days after fifth death” y de “stringency index 100 days after fifth death”

ggplot(COVID_BGC, aes(x= dpm_100d, fill = "b"))+geom_density()+
  theme_bw()+theme(legend.position = 0, panel.grid = element_blank())+
  ggtitle("Distribución de muertes por millón a los 100 días de la quinta muerte")

ggplot(COVID_BGC, aes(x= si_100d, fill = "b"))+geom_density()+
  theme_bw()+theme(legend.position = 0, panel.grid = element_blank())+
  ggtitle("Distribución del stringency index por millón a los 100 días de la quinta muerte")

Mediante el gráfico de densidad podemos ver la distribución de dos variables continuas como son el número de muertes por millón y el stringency index. Podemos comprobar la distribución de una variable binomial como es la presencia o no de campaña de vacunación en 9 periodos de 5 o 10 años.

ggplot(COVID_BGC, aes(x= vaccination_2020_2015, fill = "b"))+geom_bar()+
  theme_bw()+theme(legend.position = 0, panel.grid = element_blank())+
  ggtitle("Presencia de una campaña de vacunación entre 2015 y 2020")

ggplot(COVID_BGC, aes(x= periods_with_vaccine, fill = "b"))+geom_bar()+
  theme_bw()+theme(legend.position = 0, panel.grid = element_blank())+
  ggtitle("Número de décadas o lustros con campaña de vacunación activa")

cormat <-
    cor(COVID_BGC[2:ncol(COVID_BGC)]%>%na.omit())
cormat2 <- cormat
cormat2[upper.tri(cormat2)] <-
    NA #Para visualizar solamente una vez las correlaciones
cormat2 <- melt(round(cormat2, 2)) #Formato para poder usar ggplot
ggplot(cormat2, aes(x = Var1, y = Var2, fill = value)) + geom_tile() + scale_fill_continuous(type = "viridis")

fig <-
    plot_ly(
        x = colnames(cormat),
        y = colnames(cormat),
        z = cormat,
        type = "heatmap"
    )

fig
ggplot(COVID_BGC,
       aes(x = dpm_100d, y = vaccination_2020_2015, label = country_name)) +
    geom_jitter(position = position_jitter(seed = 1)) +
    geom_label_repel(size = 2, position = position_jitter(seed = 1)) +
    xlim(c(-100, 800))

ggplot(COVID_BGC,
       aes(x = dpm_100d, y = last_vaccine_year, label = country_name)) +
    geom_jitter(position = position_jitter(seed = 1)) +
    geom_label_repel(size = 2, position = position_jitter(seed = 1)) +
    xlim(c(-100, 800))

ggplot(COVID_BGC,
       aes(x = dpm_100d, y = first_vaccine_year, label = country_name)) +
    geom_jitter(position = position_jitter(seed = 1)) +
    geom_label_repel(size = 2, position = position_jitter(seed = 1)) +
    xlim(c(-100, 800))

model <- lm(dpm_100d ~ vaccination_2020_2015 + first_vaccine_year + last_vaccine_year, COVID_BGC)
summary(model)
## 
## Call:
## lm(formula = dpm_100d ~ vaccination_2020_2015 + first_vaccine_year + 
##     last_vaccine_year, data = COVID_BGC)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -319.96  -66.61  -17.95   12.25  385.72 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           16054.899   4441.476   3.615 0.000679 ***
## vaccination_2020_2015   -44.886     61.859  -0.726 0.471328    
## first_vaccine_year       -2.499      1.024  -2.441 0.018081 *  
## last_vaccine_year        -5.465      1.597  -3.421 0.001221 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 127.5 on 52 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.4378, Adjusted R-squared:  0.4054 
## F-statistic:  13.5 on 3 and 52 DF,  p-value: 1.241e-06